Mon. Not. R. Astron. Soc. 000, 1-16 (1994) Printed 2 February 2008 (MN WF$L style file vl.4) 



On the Life and Death of Satellite Haloes 



Giuliano Taffoni 1 , Lucio Mayer 2 , Monica Colpi 3 , Fabio Governato 2,4 

1 SISSA, via Beirut 4 - 34014 Trieste, Italy, taffoni@sissa.it 

2 Department of Astronomy, University of Washington, Seattle, WA, 98195 USA, mayer@astro.washington.edu 

3 Dipartimento di Fisica, Universita Degli Studi di Milano Bicocca, Piazza delta Scienza 3, 1-20126 Milano, Italy, colpi@mib.infn.it 
4 Osservatorio Astronomico di Brera, via Brera 28, 20121 Milano - Italy, fabio@astro.washington.edu 



Submitted to MNRAS, September 2002 



ABSTRACT 

We study the evolution of dark matter satellites orbiting inside more massive 
haloes using semi-analytical tools coupled with high-resolution N-Body simulations. 
We select initial satellite sizes, masses, orbital energies, and eccentricities as predicted 
by hierarchical models of structure formation. Both the satellite (of initial mass M s ,o) 
and the main halo (of mass Mh) are described by a Navarro, Frenk & White density 
profile with various concentrations. 

We explore the interplay between dynamical friction and tidal mass 
loss/evaporation in determining the final fate of the satellite. We provide a user- 
friendly expression for the dynamical friction timescale Tdf,ivve an d for the disruption 
time for a live (i.e. mass losing) satellite. This can be easily implemented into exist- 
ing semi-analitycal models of galaxy formation improving considerably the way they 
describe the evolution of satellites. 

Massive satellites (M S) o > O.lMh) starting from typical cosmological orbits sink 
rapidly (irrespective of the initial circularity) toward the center of the main halo where 
they merge after a time Tdf, r i g , as if they were rigid. Satellites of intermediate mass 
(O.OlMh < M Sj0 < O.lMh) suffer severe tidal mass losses as dynamical friction reduces 
their pericenter distance. In this case mass loss increases substantially their decay 
time with respect to a rigid satellite. The final fate depends on the concentration of 
the satellite, c s , relative to that of the main halo, Ch- Only in the unlikely case where 
c s /ch ~ 1 satellites are disrupted. In this mass range, Tdf,ii V c gives a measure of the 
merging time. Among the satellites whose orbits decay significantly, those that survive 
must have been moving preferentially on more circular orbits since the beginning as 
dynamical friction does not induce circularization. Lighter satellites (M s> o < O.OlMh) 
do not suffer significant orbital decay and tidal mass loss stabilizes even further the 
orbit. Their orbits should map those at the time of entrance into the main halo. 

After more than a Hubble time satellites have masses M s ~ 1 — 10%M Si o, typically, 
implying M s < O.OOlAfh for the remnants. In a Milky Way like halo, light satellites 
should be present even after several orbital times with their baryonic components 
experimenting morphological changes due to tidal stirring. 

They coexist with the remnants of more massive satellites depleted in their dark 
matter content by the tidal field, which should move preferentially on tightly bound 
orbits. 

Key words: dark matter — galaxies: kinematics and dynamics — galaxies: interac- 
tions — methods: numerical and analytical 



1 INTRODUCTION. 

In the current view, structure formation in the Universe pro- 
ceeds through a complex hierarchy of mergers between dark 
matter haloes, from the scale of dwarf galaxies up to that of 
galaxy clusters. Galaxy formation occurs within dark mat- 



ter haloes while these evolve and grow through a series of 
mergers. During the assembly of these systems, various pro- 
cesses, like morphological transformations of the stellar and 
gaseous components are expected to happen. Therefore, un- 
derstanding the dynamical evolution of dark matter haloes 
is a fundamental step of any theory of galaxy formation. 
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N-body simulation are widely used to study the dynam- 
ical evolution of cosmic structures (Governato et al. 1999) 
and they have been the most useful tool to address this 
problem, so far. Studying in detail the internal dynamical 
evolution of many haloes requires however a high number 
of particles in order to resolve substructure avoiding its ar- 
tificial evaporation (Moore, Katz & Lake 1996; Ghigna et 
al, 1998; Moore et al. 1999; Lewis et al. 2000; Jing & Suto 
2000; Fukushige & Makino 2001). The heavy computational 
burden associated with such cosmological runs limits the 
level of detail at which the evolution of the internal struc- 
ture of satellites can be followed. On the other end, non- 
cosmological simulations at very high resolution, for a lim- 
ited number of systems, have shown that merging and other 
interactions between haloes (and eventually between their 
embedded luminous galaxies) like harassment and tidal stir- 
ring can dramatically affect their global properties and their 
internal structure (Huang & Carlberg 1997; Naab, Burkert 
& Hernquist 1999; Velazquez & White 1999; Moore et al. 
1996, 1998; Mayer et al. 2001a,b, Zhang et al. 2002). 

A different approach to the problem of structure forma- 
tion and evolution is brought about by semi-analytical meth- 
ods. The backbone of the semi-analytical models of galaxy 
formation (Somerville & Primack 1999; Kauffmann et al. 
1999; Cole et al. 2000) is the merging history of dark mat- 
ter haloes which can be Monte-Carlo generated (Somerville 
& Kolatt 1999; Sheth & Lemson 1999; Cole et al. 2000) or 
calculated from N-Body simulations (Kauffmann, White & 
Guiderdoni, 1993). The evolution of substructures in semi- 
analytical models is followed in a simplified way: a merging 
event between unequal mass haloes takes place when the 
lighter halo reaches the center of the more massive one. The 
time scale for this to happen is obtained from the local ap- 
plication of Chandrasekhar's formula (1943) for dynamical 
friction. 

However, as the magnitude of the frictional drag de- 
pends on the mass of the satellite and this is a time- 
dependent quantity, we expect stripping to ultimately affect 
the orbital decay rate. Somerville & Primack (1999) include 
a simple recipe which accounts for mass stripping, reducing 
the mass of the satellite by re-calculating its tidal radius 
while it spirals toward the center along a circular orbit. 

Colpi, Mayer & Governato (1999, hereafter CMG99) 
quantified the interplay between dynamical friction and tidal 
stripping for a selected sample of orbits and satellites mass. 
Using high-resolution N-Body simulations, they showed that 
small satellites (with initial masses 50 times smaller than 
that of the primary) undergo tidal mass loss and their orbit 
decay as if they had an "effective mass" ~ 60 per cent lower 
than the initial; on typical cosmological orbits they never 
merge at the center of the primary because the magnitude 
of the drag is drastically reduced at such a small effective 
mass. The fraction of mass lost by the satellite is strictly 
related to the particular orbital parameters and halo profile 
assigned to the haloes. In order to improve this recipe and 
make it more physically motivated it is necessary to recog- 
nize that mass loss is the consequence not only of the initial 
tidal truncation but also of the repeated gravitational shocks 
occurring at each pericenter passage (Taylor & Babul 2000, 
TB; Gnedin, Hernquist & Ostriker 1999, GHO; Weinberg 
1994); the strength and effectiveness of the shocks depends 
on the central density profile and orbit of the satellite and 



might lead to its complete disruption before the merger is 
completed. This regime of disruption or, at least, of mass 
evaporation during orbital decay, is completely neglected 
by semi-analytical models of galaxy formation, even by the 
recipe adopted by Somerville & Primack (1999). However 
it has been shown that satellite orbits are very eccentric in 
CDM models (e.g. Tormen 1997; Ghigna et al. 1998), and 
this points to strong tidal shocks. 

The full dynamical evolution of the satellites must be 
studied using haloes similar to those forming in cosmolog- 
ical simulations, that have cuspy density profiles. In this 
work we will consider haloes with Navarro, Frenk & White 
(1996; 1997; hereafter NFW) profiles as opposed to the pre- 
vious work, where our analysis was restricted to isother- 
mal spheres with cores (CMG99). We note that more recent 
higher resolution simulations (Moore et al. 1999; Ghigna et 
al. 2000; Bullock et al. 2000; Jing & Suto 1999; Governato, 
Ghigna & Moore 2001) find that the inner slope of the den- 
sity profile is even steeper than the NFW. 

Following the same philosophy of CMG99, we use semi- 
analytical methods to describe the orbital evolution and 
mass loss of satellites in an NFW profile, and we compare 
the results with high resolution N-Body simulations. In par- 
ticular, we will use the theory of linear response to model dy- 
namical friction and study orbital decay (Colpi 1998; Colpi 
& Pallavicini 1998). We will apply the theory of gravitational 
shocks developed by GHO to model tidal mass loss and the 
disruption of satellites (Taylor & Babul 2001; Hayashi et al. 
2002). 

The paper is organized as follow: we first review the 
main features of NFW haloes (Section 2), and of the drag 
force as derived using the theory of linear response (Sec- 
tion 2.1). In Section 3 we study the orbital decay of a rigid 
satellite. We then move on describing the effects of the tidal 
perturbation both when the orbit is stable and when it de- 
cays due to dynamical friction (Section 4). Finally we discuss 
the global effect of dynamical friction (DF) and mass loss 
on the evolution of satellites. 



2 ORBITAL EVOLUTION IN A NFW 
PROFILE. 

A realistic representation of the density profile consistent 
with the findings of structure formation simulations is 
needed for a meaningful study of the disruption of satel- 
lites. Here, we use the so called 'universal density profile' of 
Navarro, Frenk & White (1996): 

p(r) = Jh <k (1) 

where x = r/R b is the dimensionless radius in units of the 
virial radius Rh, Mh is the mass of the halo inside Rh, Ch = 
r s /Rh is the concentration parameter (r s is a scale radius), 
and S c = Ch/[ln(l + c h ) - c h /(l + c h )]. 

A halo of given mass and size does not have a unique 
NFW profile; the concentration c plays the role of a free 
parameter that basically tells how much of the total mass is 
contained within a given inner radius. Haloes with a higher 
concentration have more mass in the central part and should 
thus be more robust against tidal effects. We will consider 



© 1994 RAS, MNRAS 000, 1-16 



On Life and Death 3 



various concentrations for both the primary halo and the 
satellite. 

The mass profile of a spherically symmetric halo (i.e. the 
mass contained inside a sphere of radius r) can be obtained 
integrating equation (1) over the spherical volume 



M(r) = M h 



ln(l + Chi)- c h x/(l + c h x) 
ln(l + c h ) -Ch/(1 + Ch) 



(2) 



and used to calculate the circular velocity profile, V c 2 (r) = 
GM(r)/r, and the one-dimensional velocity dispersion a(r) 



a 2 (r) = 75.53 V c 2 (2.157? h /c h ) (c h x)(l + c h x) 2 l(c h x) 



(3) 



X(x) 



ln(l + y) 1 

y 3 (l+y) 2 y 2 (l + y) 3 



dy 



(Kolatt et al. 2000). 

The gravitational potential of a NFW halo can be writ- 
ten as: 

m = -K? ( r) + Y!?:ttt ft f . w 

ln(l + c h ) - c h /(l + c h ) 

here Vh is the value of the circular velocity at the virial radius 
Rh. The orbits in this potential can be determined using 
the planar polar coordinates r(t) and 9(t), solving for the 
equation of motion (Binney & Tremaine 1987). The motion 
of a satellite is then determined by the initial specific angular 
momentum J and orbital energy E, or equivalently by the 
radius r c (E) of the circular orbit having the same energy E, 
and by the circularity e = J/J c , where J c = Vh(r c ) • r c (E). 
We define a generalized orbital eccentricity: 

(5) 



po 7 pc 



here r apo and r pcr are the roots of the orbit equation: 



1 

— + 

r 2 



~ E] 

J 2 







(6) 



(Binney & Tremaine 1987), and they are respectively the 
apocenter and the pericenter radii of the orbit. Using the 
previous equation it is possible to derive a relation between 
e and the orbital parameters , so that for each value of 
r c (E) and e we can determine the apocenter and pericenter 
distances of the orbit. We introduce the dimensionless radius 
of the circular orbit x c (E) = r c (E)/Rh- 

From equation (6), the orbital period is 



Porb = 2 



dx 



^2[0.5K 2 (r c ) + <p(rc) - tfx)] - J 2 / 'a 



(7) 



A satellite, described by a NFW profile, has subscript s in 
all the corresponding halo properties. The satellite, before 
mass loss, has a mass M Sj o inside its virial radius -Rs.o- 
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Figure 1. Dynamical friction timescale Tdf jr i g versus M B fi/M^ 
for a satellite in a Milky Way like halo with x c (E) = 0.5 and 
e = 1; filled symbols are from TLR, while open symbols are from 
the local approximation of dynamical friction as given by solving 
equation (9). Dots refer to Ch = 7 while squares for Ch = 14; the 
solid line corresponds to the fit Tdf, r ig ~ 1.3R^VhX 2 /[GM s fi In A] 
where InA = ln(l + M h /M B>0 ). 



Weinberg 1989 for a study of dynamical friction in a self- 
gravitating medium). The dissipative force on the satellite 
is computed tracing, in a self-consistent way, the collective, 
global response of the particles to the gravitational pertur- 
bation excited by the satellite. The force includes the tidal 
deformation in the density field (absent in an infinite uni- 
form medium), the trailing density wave which is evolving 
in time, and the shift of the bary center of the primary. We 
omit here the complex expression of the force referring to 
Colpi and Pallavicini (1998) and CMG99 for details. We 
only remark that the frictional force at the current satellite 
position R(t) can be written formally as an integral upon 
space and time 



F DF (t) = —GM t 



dt'Apftt-t')-^ 



(8) 



\R(t)-m\ 3 

where Ap(r, t — t') maps of the time dependent disturbances 
in the density field created, over time, by the satellite in 
its motion. We estimate the force specifically for the NFW 
density profile. TLR do not contain any free parameter (no 
Coulomb logarithm) except the mass and the radius of the 
satellite, equation (8) describes the sinking of satellites mov- 
ing along orbits of arbitrary eccentricity even outside the 
primary halo. 



2.1 The Theory of Linear Response 

The theory of linear response (hereafter TLR) is a relatively 
novel approach to the study of dynamical friction in the non- 
uniform stellar background of a spherical self-gravitating 
halo (Colpi & Pallavicini 1998; Colpi 1998; CMG99; see also 

* For an Singular Isothermal Profile (SIP) the eccentricity is only 
a function of e (van den Bosch et al. 1999) 



3 NUMERICAL SIMULATIONS 

The simulations, whose results are analyzed in the forth- 
coming sections of this paper, have been performed with 
PKDGRAV, a fast parallel binary treecode widely used to 
study structure formation and galactic dynamics (e.g. Power 
et al. 2002; Ghigna et al. 1998; Mayer et al. 2001a,b; Dika- 
iakos & Stadel 1996; Stadel, 2002). The force calculation is 
performed using a binary tree and using the Barnes-Hut cri- 
terion for evaluating the multipoles up to the hexadecapole 



© 1994 RAS, MNRAS 000, 1-16 



4 G. Taffoni et al. 



O 4 



t- 2 




M/M h =0.1 



I I I I I I I I I I I I I I I I I I I I — L 



0.2 0.4 0.6 0.8 1 



Figure 2. Dynamical friction time Tdf^ig versus circularity e for 
for M s /M h = 0.02, 0.05 and 0.1, here x c (E) = 0.5. Points arc the 
TLR data and the solid lines are the model results Tdf jr i g oc e° 
where Tdf jr i g is given by equation (12) and we estimate a using 
equation (15). 



order; an opening angle 6 = 0.7 was used in all the runs. 
The code uses a leapfrog integrator and has multistepping 
capabilities. The simulations employing a rigid satelllite use 
an NFW primary halo resolved by either 100.000 or 1 mil- 
lion particles; the satellite is modeled as in Van den Bosch et 
al. (1998) and CMG99, i.e it is represented by a point mass 
softened using a spline kernel (the same kernel adopted for 
all the particles in the simulations). The softening of the 
particles in the primary halo is 200 pc; that of the rigid 
satellite is 3.5 kpc in the reference case where the latter has 
a mass M Si o = 0.05Mh. Satellites of different masses have a 
softening scaled ~ M^q 3 . In the simulations for "live" (i.e., 
mass losing) satellites, these are resolved by either 20.000 
or 50.000 particles (the same resolutions holds in runs of 
live satellites moving in a fixed external potential) while the 
primary halo is resolved using 100.000 particles. The par- 
ticle softening for satellites of different masses is rescaled 
as for the rigid satellites and the same scaling holds also 
between these particles and those of the primary halo (as 
a reference, for a satellite of 0.05 Mh the softening is 74.4 
pc). We note that the softening of the rigid satellite is fixed 
in such a way that a deformable satellite having high con- 
centration (c = 20) has roughly the same half mass radius 
than the corresponding rigid satellite of the same mass; this 
ensures that a comparison between runs with rigid and de- 
formable satellites is meaningful (the decay rate depends on 
the softening of the rigid satellite at a given mass; see e.g. 
Van Albada 1987). By comparing runs with different resolu- 
tions for the primary halo we verified the robustness of the 
obtained orbital decay rate in absence stripping. By com- 
paring runs having deformable satellites with different reso- 
lutions we tested whether artificial heating due to two-body 
collisions was playing a role in the determining the actual 
mass loss rate. The simulations employed timesteps as small 
as 10 5 years in the inner regions of the halos, namely more 
than an order of magnitude smaller than the local orbital 



time; as a result, energy was conserved to better than 1 per 
cent. 



4 THE SINKING OF A RIGID SATELLITE IN 
A NFW PROFILE 

In this Section we explore the evolution of a rigid satellite of 
mass M Sj o orbiting inside a halo with NFW density profile, 
using TLR. The halo is scaled to the Milky Way mass Mh = 
10 12 Mq, has a tidal radius Rh — 200 kpc and concentration 
Ch = 7 or 14, within the spread of cosmological values (Eke, 
Navarro & Steinmetz 2001). 

Fig. (1) shows the dynamical friction time Tdf,ri g as a 
function of the satellite mass M Sj o (expressed in units of 
Mh); the satellite moves on a circular orbit at x c (E) = 0.5 
We find no significant dependence of Tdf, r i g on the halo con- 
centration. The fit in the figure tries to single out the depen- 
dences of r d f, r ig on the mass of the satellite and its initial 
orbit in a simple way and ties to the familiar expression 
of the dynamical friction timescale, derived in the local ap- 
proximation, for the case of an isothermal sphere (Binney 
& Tremaine 1987). If one again uses the expression of the 
frictional force, given by Chandrasekhar (1943), treating the 
background density and dispersion velocity as local quanti- 
ties (evaluated at the satellite current position) ' , the evolu- 
tion equation of a satellite spiraling down on circular orbits 
in a NFW main halo is 



l d[rV c (r)] 
r dt 



= -47rlnAG 2 M s , 



P(r,c h ) 
V c 2 (r) 



erf(Y) 



(9) 



where Y — V c (r)/\/2o"(r). This equation can be integrated 
grouping all quantities depending on r, on the left hand side 
of equation (9) to give 



/' 



®{x, Ch) dx 



GM Si o In A 



Tdf.i 



(10) 



where x c the initial radius of the circular orbit. 

The function @(x,Ch) has an analytical expression that 
can be fitted, with an average error of one part over 1000, 

as 

e(^, C h)~/( C h)^ - 97 , (ii) 

leading to a dynamical friction timescale for circular orbits 



Tdf,rig ~ 0.6/(0.) ■ 



for e = 1 , 



GM S , In A 
where In A = ln(l + M h /M Si0 ) and /(c h ) is 
/(c h ) = 1.6765 +0.0446 c h . 



(12) 



(13) 



This simple analysis explains while a fit similar to that 
for a singular isothermal sphere (see Fig. (1) and its caption) 
is acceptable even in a NFW profile. 

t Dynamical friction is the result of a "long range" disturbance 
(Hernquist & Weinberg 1989; Colpi & Pallavicini 1992; CMG99). 
Treating it as local is conceptually incorrect. The error that re- 
sults, which is difficult to quantify unless the whole treatment is 
included, is customarily absorbed in the Coulomb logarithm. 
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Figure 3. The residual mass of a satellite at the first pericenter 
passage as function of the orbital parameters, when c s /ch = 2. We 
assume that the "tidal cut" instantaneously reduces the satellite 
mass: Each curve is labeled with the residual mass in units of the 
initial one. The vertical dotted line is the most probable value 
of the eccentricity in a cosmological environment (Tormen 1997), 
the dashed vertical lines are the lcr variance. 



The frictional time is also a function of the initial orbital 
circularity e: We thus explored the dependence of Tdf,rig(e) 
as accretion of satellites in a main halo occurs preferentially 
along rather eccentric orbits. This dependence has been al- 
ready discussed in Lacey & Cole (1993), van den Bosch et 
al. (1999), and CMG99 for isothermal profiles, giving 

T"df,ri g (e) ~ T d f,rig(e = l)e a . (14) 

CMG99 noticed that, for a fixed satellite mass (M B ,o <C Mh), 
the timescale varies with the orbital energy, suggesting for a 
a dependence on x c (E) (see CMG99 for the suggested values 
of a). 

In a NFW halo, we find that a depends on x c {E), and on 
M s fi/Mh- We find that, whereby relatively heavy satellites 
decay on a time almost independent of e, lighter satellites 
decay on much shorter times when e — > 0. This is shown in 
Fig. (2). A useful fit to a as a function of orbital energy and 
mass ratio is 

a{x c ,M sfi /M h ) ~ 
0.475 <! 1 - tanh 



10.3 



7.5x c 



(15) 



As in CMG99, we used high resolution N-body simulations 
to follow the orbital decay of a rigid satellite (using up to 10 6 
particles to sample the halo mass distribution) and found a 
very close match with the theory of linear response. The 
two approaches agree in a number of details on the evolu- 
tion, the most remarkable being the temporary rise of the 
orbital angular momentum observed during the final stages 
of the decay. This is a manifestation of the fact that in the 
background medium, no longer uniform, the satellite moves 
inside or close to its distorted wake that, near pericenter 
distance, induces a positive torque (Colpi et al. in prepara- 
tion). 

As a final remark we notice that initial eccentric orbits 
that decay by dynamical friction do not chance significantly 



their degree of circularity with time, as was also found in 
isothermal profiles (van den Bosh et al. 1999; CMG99). 



5 THE DYNAMICAL EVOLUTION OF A LIVE 
SATELLITE 

The evolution of a rigid object is determined by the frictional 
drag force and its survival time corresponds to the dynami- 
cal friction time Tdf.rig- However a real satellite is not a rigid 
point mass but a deformable distribution of particles mov- 
ing inside a halo. Its life is then dramatically influenced by 
the tidal perturbations induced by the gravitational field of 
the primary halo. The global effect of the tidal perturbation 
is the progressive evaporation of the satellite. This process 
takes place during the orbital evolution and it is generally 
sensitive to the internal properties of the satellite and of the 
surrounding halo. 

Our aim is to model realistically the tidal effects in order 
to evaluate the mass that remains bound to the satellite, 
M s (t), each time along the orbit. 

We distinguish two tidal effects: a tidal truncation (tidal 
cut), originated by the average tidal force exerted by the 
main halo at the distance of the satellite, and an evapora- 
tion effect induced by the rapidly varying tidal force near 
pericenter radii for satellites moving on an eccentric orbit. 
In the latter case we speak of tidal shocks - short impulses 
are imparted to bound particles within the satellite, heating 
the system and causing its dissolution. 

5.1 The tidal truncation 

A tidally limited satellite is truncated at its tidal radius 
-Rs.tid, which, loosely speaking, corresponds to the distance 
(relative to the satellite center) at which the mean density 
of the satellite is of the order of the mean density of the 
hosting halo, at the satellite position r: 



(16) 



The evaluation of the tidal radius requires a relation 
between 7? Sit id and r which is customarily derived from the 
force equivalence between internal gravity and external tides 
leading to the implicit equation: 



M S (R S 



(2 - d In M h /d In r) M h (r) 



1/3 



(17) 



(Tormen, Diaferio & Syer 1998). The mass tidally lost, 
AM Sjt id is thus computed subtracting spherical shell above 
^s.tid, using equation (2) . While strictly valid for a satel- 
lite moving on a circular orbit (where the combined potential 
over the system is static in the satellite frame) i? ajt id gives, if 
evaluated at every single point r (Binney & Tremain 1987), 
an approximate expression for the instantaneous tidal ra- 
dius in the case of non circular motion. This implies that, 
on stable orbits, AM Sj tid is maximum at the first pericenter 
passage; the mass of the satellite would then remain con- 
stant. In Fig. (3) we give the residual mass after instanta- 
neous tidal cut, as a function of circularity, as computed 
using equation (17). 

Tidal stripping however does not occur instantaneously, 
and, following TB suggestion, we model mass loss, over a few 
orbital periods, adopting the expression 
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dM 
dt 



AM s , tid (t) 
2tt/u;(£) 



(18) 



where u(i) is the instantaneous orbital angular velocity. This 
is compared with results from numerical simulations. Fig. (5) 
gives the satellite mass as a function of time for a selected 
run. We find that mass loss by tidal cut, as described by 
equation (18), reproduces the result of our N-body simula- 
tion only in the early phase: the satellite loses mass at a 
peace larger than predicted by equation (18) (we refer to 
the dashed line of Fig. [5]). We believe that this is due to 
the action of tidal shocks (and not a numerical artifact). 

The number of particles in the N-body simulations is 
chosen in order to avoid as much as possible numerical two- 
body relaxation which could increase, artificially, the overall 
evaporation rate (Gnedin & Ostriker 1999; Moore, Katz & 
Lake 1996) . Numerical relaxation disperses satellite parti- 
cles over a timescale related to the number of particles N 



trh = 0.138 



G°- 5 m 4 ln(0.4A0 



(19) 



where 7?hm is the half mass radius and m» is the particle 
mass, and M s = m*N. As shown in Table (1) the initial 
relaxation time is ~ 100 Gyr and remains longer than 10 Gyr 
as mass loss continues. This is an indication that numerical 
two-body relaxation is unimportant. We thus proceed on 
modeling mass loss with the inclusion of tidal shocks. In 
the next Section we estimate the shock-induced evaporation 
time and show that, for NFW satellites, it can be shorter 
than cosmic age. 



5.2 Heating &z evaporation 

The description of the dynamical evolution of a satellite 
must include also tidal heating due to compressive tidal 
shocks. 

The theory of shock heating was developed by Ostriker 
ct al. (1972) and Spitzer (1987) to model the evolution of 
globular clusters. Recent works by Gnedin & Ostriker (1997) 
and GHO extend this theory also to tidal perturbation on 
satellites moving on eccentric orbits inside an extended mass 
distribution. We will use the GHO model to treat tidal 
shocks on satellites orbiting specifically inside a NFW halo. 
At each pericenter passage satellites cross the denser regions 
of the main halo: the rapidly varying tidal force induces a 
gravitational shock inside the satellite. The shock increases 
the velocities of satellite particles, and reduces the satellite 
binding energy. As a result the satellite expands. 

The amount of heating is a function of the orbital pa- 
rameters and of the concentration of the main halo: 



(AS) = T (c h ,x c [E],e) A(x T ) ■ R B 



(20) 



The shock is more intense in the outer layers, as it depends 
also on the satellite radius R B (see Appendix A for the de- 
tails of the calculation). As suggested by GHO, we use an 
adiabatic correction A(x T ) = (1 + x 2 ^) 1 with 7 = —5/2 (see 
also Weinberg 1995) . Here x T = lot is the adiabatic pa- 
rameter, t is the duration of the shock and to = a a (R B )/R B 
where <r 3 is the velocity dispersion of particles in the satel- 
lite at radius R B . The value of r is related to the pericenter 
crossing time: we assume r = 0.5/w por , where u pcT is the or- 
bital angular velocity at pericenter distance. The adiabatic 
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Figure 4. The amount of shock heating as function of the circular 
radius x c (E). (AE) is normalized to the value of (AE) when 
x c (E) = 0.5 and e = 0.6. We consider three different values of 
the circularity: e = 0.4 (solid line), e = 0.6 (dotted line) and 
e = 0.8 (dashed line). The satellite and halo concentration are 
chosen such as c B /c^ = 2. 



correction accounts for the fact that the susceptibility of a 
system to the tidal shocks will also depend on its internal 
dynamics - when the internal orbital time is very short a 
particle in the satellite will receive two opposite tidal kicks 
of nearly the same magnitude and the net effect will be small 
(e.g. Binney & Tremaine 1987). 

We introduce a characteristic shock timescale com- 
puted, after each pericenter passage, as 



t B h' — 



Porb Eq 

~ ' (AE hm ) ' 



(21) 



where Eo — 0.25GM SjPer //? Si hm is the binding energy of the 
tidally truncated satellite of mass M BiPCT evaluated accord- 
ing to equation (18) at the time of pericentric passage. Both 
Eo and (A£hm) are evaluated at the half mass radius R B ,hm 
which is a function of the satellite concentration. A second 
order energy change due to shock heating is responsible for 
increasing the internal velocity dispersion and allows ad- 
ditional particles to leave the satellite. To account of this 
second order perturbation we assume that i s h = 0.43t s h' 
(Gnedin, Lee & Ostriker 1999; hereafter GLO). Table (1) 
shows the shock time for the satellite modeled at first peri- 
center passage. The number of pericenter passages roughly 
necessary to unbind the satellite is t B h/Poib- Lastly we notice 
that (AE) increases linearly with the halo concentration Ch, 
because in highly concentrated haloes the gradient of the 
gravitational force is steeper. 

The amount of heating is also a function of the or- 
bital parameters; in Fig. (4) we study the energy gain as a 
function of x c (E) for different values of the circularity. The 
fast growth of {AE) for small values of x c (E) confirms that 
shocks on radial orbits are more intense: a satellite moving 
on a circular orbit is not subject to any heating. 
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Table 1. The characteristic timescalcs 

Model c s /c h t sh [Gyr] P orb [Gyr] t rh [Gyr] 

Low concentration 

e = 0.7 x c = 0.5 0.5 12.6 4.7 176.4 

e = 0.5 x c = 0.3 0.5 0.7 2.6 173.6 
Intermediate concentration 

e = 0.7 x c = 0.5 1 93.6 4.7 119.2 

e = 0.5 x c = 0.3 1 2.00 2.6 112.6 
High concentration 

e = 0.7 x c = 0.5 2 130 4.7 80.7 

e = 0.5 x c = 0.3 2 6.6 2.6 73.8 
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Figure 5. Bound mass in units of the initial mass as a function 
of time, for a satellite moving on a stable orbit. The orbital pa- 
rameters are e = 0.65 and x c (E) = 0.5. The concentration ratio 
is c s /ch = 0.7 (top panel) and c B /c^ = 0.4 (bottom panel). The 
symbols are the N-body data and the solid line the semi-analytical 
model. Stars identify each pericenter passage. The dashed line is 
the bound mass that would remain if we apply only the tidal cut 
using equation (18) over a few orbital periods. 



5.3 Modeling the mass loss 

Tidal shocks are events leading to the escape of particles. To 
model the induced mass loss, we introduce the so called es- 
cape probability function £ s h, analog of £ c customarily used 
to describe globular cluster evaporation by two-body relax- 
ation processes (Spitzer 1987). 

Mass loss can be predicted using the dimensionless rate 
of escape 

M(t) dt ■ 
Similarly here we define 
U dM 

^ sh ~ M(t) dt ' (26) 

For the case of escape by two-body relaxation £ s h is a con- 
stant known to vary from 7.4 • 10~ 3 for an isolated halo to 
0.045 for a tidally truncated halo (Spitzer 1987). On the 
contrary, when tidal shocks are present and dominate, the 
escape probability becomes a function of time; £ s h peaks just 
after each pericenter passage (GLO) , rapidly decreasing un- 
til the next shock event. £ a h is then a periodic function of 
period P orb and we find that it can be fitted using both 
simulations and results by GLH as: 

UWoc(^)- a5 exp-(^)°' 5 , (24) 

where t tT — 13 i s h and t pcl is the pericenter time. The shock 
escape probability is equal to unity at t = t pcr + t dyn where 
tdyn is the dynamical time of the satellite. The shock time 
must be evaluated at each pericenter passage as it varies 
according to the current orbit and mass of the satellite. 

The orbital timescale is sometime shorter that the shock 
timescale so that the satellite suffers. If it becomes unbound, 
further dispersal of the last particles occurs on the crossing 
time of the damaged system. 



5.4 Testing the model for a live satellite 

The dynamical evolution of a satellite is described using 
a semi-analytical code which accounts of both dynamical 
friction and mass loss. In this context, we use the expression 
of the drag force as given in equation (9), since it is much 
faster, and in close match with TLR (see Section 3). At 
each time step we upgrade the satellite mass according to 
equation (18) and equation (23). To test the ability of our 
code to follow the evolution of an NFW satellite, we compare 
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the results with those derived from a selected set of N-body 
simulations. 

5.4-.1 Tidal perturbations on stable orbits 

To isolate and study the effect of a pure tidal perturbation 
we explore the dynamical evolution of a low concentration 
satellite on an unperturbed orbit. In this case, the heating by 
tidal shocks varies solely as a consequence of the progressive 
reduction of the satellite half-mass radius. For this reason 
we expect a progressive reduction of the shock destructive 
power as time passes. In Fig. (5) we show the evolution of a 
satellite with mass M Sj o = O.OlMh; the orbital parameters 
are chosen to reflect a typical cosmological orbit: e = 0.65 
and Xc{E) = 0.5. The bottom panel shows a low c s /c h satel- 
lite, disrupted after the second passage to pericenter. The 
top panel shows the evolution of a higher c s /ch satellite sur- 
viving for more than 12 Gyr, despite having lost more than 
the 99 per cent of its mass. Until the first pericenter pas- 
sage only the tidal cut accounts for the mass loss. The good 
agreement between the simulation and the code before the 
first pericenter passage suggests that the TB recipe is accu- 
rate enough to reproduce the mass loss before (or in absence 
of) the shock heating. 
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Figure 7. We have plotted the ratio of the dynamical friction 
time of a non-deformable satellites to the same time for a live 
satellite of equal mass, initially, as function of M^o/M^. The 
concentration ratio is c B /ch = 2, and x c (E) = 0.5. We vary the 
circularity which is labeled with different symbols. We notice that 
symbols for e = 0.3 refer to a different orbital energy: x c (E) = 0.3. 



5.4-2 The combined effect of dynamical friction and tidal 
stripping 

The dynamical evolution of a satellite is driven by the com- 
bined effect of dynamical friction that drives the satellite 
to the center of the main halo, and the tidal perturbation 
which reduces its mass. The two processes are intimately 
connected as the drag force is strongly related to the mass 
and size of the satellite. 

In Fig. (6) we compare the semi-analytical model with 
the the results of N-body simulations for satellites with 
Cs/ch = 2. The initial orbital parameters are e = 0.7, and 
x c (E) = 0.5. We study two different cases: a light satellite of 
mass M Sj o = 0.02 Mh and a massive one with M B = 0.1 Mh. 
The mass loss rate and the orbital evolution are well repro- 
duced in both cases. The massive satellite loses mass during 
evolution, yet a core of bound particles survives having 5% of 
its initial mass but sinks to the center merging with the main 
halo in 3 orbital periods. On the contrary, the light satellite 
loses 99% of its mass but a bound core remains which moves 
on an inner orbit stable against dynamical friction, following 
mass loss. 



6 THE FATE OF SATELLITES 
6.1 The recipe 

We now use our semi-analytical model to quantify how mass 
loss affects the orbital decay. In Fig. (7) we give as a func- 
tion of Msfi/Mh, the ratio of the dynamical friction time of 
a rigid satellite Tdf, r i g to the same time for a homologous live 
satellite Tdf,n V e ■ In taking the ratio we mainly quantify the 

t We use highly concentrated satellites (c s /ch = 2) that are not 
rapidly disrupted by tidal interactions. 



importance of mass loss in affecting the lifetime of a satellite. 
Fig. (7) shows that massive satellites (M s ,o/Mh > 0.1) sink 
to the center of the main halo on a timescale Tdfjifc ~ Tdf.rig 
as if they were rigid. In the mass range M Sj o — 0.005 — O.lMh 
the satellites sink toward the center of the main halo by dy- 
namical friction and so lose mass efficiently. Accordingly, 
the dynamical friction time Tdf,H V e is a few to several times 
longer than r d f >rig . At lower masses, M Sj o ~< 0.005 M h , or- 
bits are less perturbed by friction and mass stripping be- 
comes less important. Thus the ratio Tdf,rig/Tdf,nfe starts to 
rise again. 

We now estimate the dynamical friction timescale for a 
live satellite in three different regimes. 

For massive satellites, M Sj o > 0.1 Mh, the dynamical friction 
time is not affected by the mass loss so 



7"df,livc ~ 7~df,rig 

where: 



°- 5 GM7 Aig 



'-,Ch,x c (E) 



A 



c h ,x c (E) = f(ch 



My 



xl- 97 (E) 



In (1 + M h /M s ,o) ' 



(25) 



(26) 



and f(ch) is given by equation (13). In this case, the dynami- 
cal friction time depends weakly on e; the exponent a ~ 0, as 
indicated by equation (15). For 0.007 M h < M Sj0 < 0.08 M b , 
we provide a fit of the form: 



Tdf.live 



Rj Vh 

GM S , 



-Aiiv 



M sfi c s 

-T7-,— ,x c (E),e 
Mh c h 



(27) 



For e = 1 



/ Mg£ C. \ _ 

" 4l,v ° V M h 'c h ' Xc ' e ~ ) ~ 



ll ' 2r > 0.1)7- • 1.123 



(c s /c h ) 6 



Ch 



(28) 



where 
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Figure 6. We plot the time evolution of the mass, in units of the initial mass (left panels), and orbital angular momentum scaled to the 
initial one (right panels). The adopted values are: M a ,o = 0.1 Mj, and M s ,o = 0.02 M^. The concentration ratio and orbital parameters 
are c s /ch = 2, e = 0.7 and x c (E) = 0.5, respectively. Points are from N-body runs and solid line from our semi-analytical model. 



B{x ) = -0.050 + 0.335 x c + 0.328 x\ , 



C(x c 



2.151 - 14.176 x c + 27.383 



(29) 
(30) 



This fit reproduces the semi-analytic estimate of the decay 
time of a live satellite on circular orbits with an error of 9 
per cent, for 0.2 < x c (E) < 1. 

For eccentric orbits we find that: 



An- 



( Ms£ 



Ch 



M h ' c h ' 



Xc, e = 1 



where 
M 



Q 



(M s , \ 



°- 4+2 (tf'-) X(e -°- 2) 



0.9 + 10 s (12.84 + 3.04x c - 23.4xc 



(31) 



/ M s ,p O0i 

V M h 1-1 



0.0 077 

08x c 



0.0362 



(32) 



This formula holds when 0.3 < x c < 0.9 and 0.3 < e < 0.8 
and reproduces the semi-analytical data within and error of 
15%. 

Interestingly, we notice that for an eccentric orbit the 



decay time can be longer that the DF time on the circular 
orbit with the same initial orbital energy, since mass loss on 
eccentric orbit is higher because of the tidal shock. 

For e > 0.8, .Alive we use the interpolation of equations 
(31) and (28). If 0.08 M h < M 8 , < 0.1 M h , we suggest to 
linear interpolate equations (25) and (27). 

Satellites M Si o 0.007 Mh, evolve on slightly perturbed 
orbits, the dynamical friction timescale in this case is al least 
two times longer than for the rigid satellite (CMG99 simi- 
larly found an increase of a factor e for M s ,o = 0.02Mh). We 
suggest to use equation (12) for Tdf.iafej in this mass range, 
increased by a factor ~ 2.73, together with equation (15) 
for the exponent a, to account for the correction to circular- 
ity. For very light satellites (M s ,o < 10 _4 M h ) equation (14) 
holds. 

In Appendix B we give a simple expression for the dis- 
ruption time tdis that can be used for a comparison with the 
other timescale. 
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Figure 9. The life diagram of a satellite with M Si o = 0.01 M\y. Each plot is labeled with the value of the x c (E). We identify the region 
of the parameters space where the satellite sinks to the center of the main halo (M), evaporates in the background (D) or survives (S). 




Figure 8. The dynamical friction time of life satellites As func- 
tion of the satellite initial mass. Symbols are the results obtained 
from the semi-analytic code, lines are the prescription of the fit- 
ting formulae (from equation [27] to equation [32]). We set the 
circularity e = 0.7 and the relative concentration c B /c^ = 2 and 
we vary the orbital energy: x c (E) = 0.5 (solid line and squares), 
x c (E) = 0.7 (dotted line and stars), x c (E) = 0.3 (dashed line and 
circles). 



6.2 Merging, distruption or survival 

We here investigate the dynamical evolution of a population 
of satellites, in a given halo. An individual satellite is labeled 
by four parameters; x c (E) and e identify the orbit, while 
initial mass M s ,o and concentration c s identify the internal 
properties. Each combination of the four parameters leads to 
a different final state for the satellite: rapid merging toward 
the center of the main halo (M), disruption (D), or survival 
(S) (when a residual mass M s remains bound and maintains 
its identity, orbiting in the main halo for a time longer than 
the Hubble time). 

An important role is played by the concentration ratio 
as shown by the life diagrams in Fig. (9). These predict the 
final fate of a satellite with M s ,o = 0.01 Mh, as function 



Figure 10. Probability distribution for the three final endpoints: 
merging (M), disruption (D) or survival (S) as function of the 
initial satellite mass. The thick solid line refers to the case of a 
rigid satellite. 



of the Cs/ch and of the orbital parameters. The fractional 
area in this parameter space leading to disruption, survival 
or decay is an estimate of the relative importance of these 
processes in determining the satellite fate. Disruption due 
to the tidal perturbation is the fate of those satellites that 
initially move on close orbits despite c s /ch- Satellites moving 
along typical (plunging) cosmological orbits survive over a 
Hubble time only if they had a concentration higher than 
that of the main halo at the time of their infall. 

In Fig. (10) we have drawn the probability distribu- 
tion relative to the three final states: direct merging (by dy- 
namical friction), which dominates at large masses, survival 
and/or disruption which is the most likely end for satellite 
with M s ,o < 0.01 Mh. In producing Fig. (10) we have gen- 
erated evolutionary paths (ending after a time equal to the 
Hubble time) for satellites starting from a flat distribution 
of orbital parameters and concentrations. 

Our study suggests that those satellites that survive 
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Figure 11. The mass that remains as a function of the orbital 
parameters after more than a Hubble time, for a satellite with 
A^s,0 = 0.01 Mh. Dots identify the regions where the satellite is 
disrupted. The contour line on the xy plane identify the loci where 
M s /M s , =0.1. 
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Figure 12. The distribution of the final mass of a satellite of 
initial mass Ms = O.OlMh after more than a Hubble time. His- 
tograms are derived starting from a uniform distribution of orbital 
parameters, for two values of the concentration ratio: c s /ch = 2 
(filled grey area) and c s /ch = 1. 



have lost memory of their initial state: dynamical friction 
perturbs the orbit and tidal stripping reduces the satellite 
mass. In Fig. (11) we compute the mass of the satellites 
that remains after a Hubble time. The figure refers to a 
high concentration case, but we extend our analysis also to 
low concentrated satellites as shown in Fig. (12), where we 
compute the cumulative mass distribution for all the initial 
orbital parameters. On average, much less then 10% of the 
initial mass remains bound. Of course, in general circular 
orbits do not cause serious damages to the satellite as shock 
heating is less intense (an exception is represented by satel- 
lites on very tightly bound orbits). In Fig. (11), dots show 



0.01 




Figure 13. The apocenters radius in units of its initial value for 
a satellite M s< q = O.OlMh evaluated after more than a Hubble 
time (we terminate at 15 Gyr). The two plots refer to c s /c h = 
2 (bottom) and c s /c^ = 1 (top). Dots identify satellites that 
evaporate before 15 Gyr. The contour lines on the xy plane select 
the regions where the relative reduction of r apo is 0.1 (dotted 
lines) and 0.5 (dashed lines). 



the final mass just prior evaporation. As expected, radial 
orbits can more easily dissolve a satellite. 

The strength of the orbital decay can be estimated mea- 
suring the reduction of the apocenter distance. In Fig. (13) 
we plot the distribution of apocenter radii for a satellite 
with M 8 ,o = O.OlAfh after 15 Gyr of orbital evolution. The 
strength of the drag force reduces the apocenter distance by 
a factor of two for cosmological orbits and it is not signifi- 
cantly affected by the concentration. 



6.3 Cosmological examples 

Now, we apply our analysis to some cosmological relevant ex- 
amples. We discuss the evolution of different satellites which 
orbit in a cluster-like and galaxy-like haloes. The cluster halo 
is a Coma-like cluster with mass M h = 5 - 1O 15 M and con- 
centration Ch = 3.44. The Milky Way-like galaxy halo has 
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M h = 1O 12 M and c h = 10.44. For all cases, the initial 
orbital parameters are chosen as e = 0.6 and x c (E) = 0.5. 

Group in Coma 

We consider a massive group-like satellite of mass M s ,o = 
3 • 10 13 Mq and c s = 7.5 which enters the Coma-like halo 
at z = 0.5. In a ACDM Cosmology it evolves for ~4.8 Gyr 
inside the halo until z = 0. 

As suggested by the high value of c s /ch the satellite 
is not disrupted. Since M s ,o = 0.006Mh the orbit is stable 
and, with this choice of the initial orbital parameters, the 
satellite evolves for ~ 1.5 P or b- The final apocenter radius 
is r apo ~ 0.85r apo ,o and its final mass is M B — 7.2 • 10 12 Mq. 

Milky Way in Coma 

A Milky Way-like satellite has mass M s , = 10 12 M© and 
c s = 10.44. If it enters the Coma-like halo at z = 0.5 it 
evolves for ~ 1.5 P or b- The orbit remains almost unper- 
turbed (r apo ^ 0.99r apo ,o) as strength of the drag force is ex- 
tremely weak since M s ,o = 0.0002Mh. Due to the extremely 
high relative concentration, c s /ch ~ 3 the satellite does not 
evaporate and its final mass ad z — is M s = 2.5 ■ lO n M0 

Large Magellanic Cloud in Milky Way 

A Large Magellanic Cloud halo has M Sj o = 10 11 M Q and 
c s = 11.9. As expected due to its relative high mass, the 
satellite merges with the Milky Way in ~ 4 Gyr. Before 
merging, the satellite loses 97% of its mass that is dispersed 
in the Milky Way Halo. 

Dwarf in Milky Way 

We consider a Dwarf-like satellite of mass M s , = 5- 1O 9 M 
and concentration c a = 13.6. If it enters the Milky Way-like 
halo at z = 0.5 it evolves on an almost unperturbed orbit 
for ~ 2 Porb and its final mass is M s = 2 • 10 s Mq. 

Dwarf in Milky Way at high redshift A Dwarf-like 
satellite enters a Milky Way-like halo at z = 2, when 
the Milky Way has mass Mh = 1O 11 M0 and concentra- 
tion c h = 6.15. The satellite has M s , = 5 • 1O 9 M and 
c s = 6.8. The dwarf evolves for ~11 Gyr. Due to its low 
relative concentration it loses 99% of its initial mass dur- 
ing the first orbital period; its orbit then becomes stable 
(r apo — 0.36r apo ,o). Note that here we are not accounting 
for the evolution of the main halo which grows in mass dur- 
ing the remaining 11 Gyr before z = 0, but the influence 
of accreted mass on the central regions dynamics should be 
relatively small (Helmi et al 2002) Finally we notice that 
our Milky Way model does not account for the presence 
of a disk. Penarrubia, Kroupa & Boily (2002) suggest that 
orbital evolution changes when the potential well has a flat- 
tened component, and DF is more efficient for satellites with 
low-inclination orbits (respect to the disk): the orbital de- 
cay is accelerated and the orbital plane decays over the disk 
plane. They also find instead that DF enhances the survival 
time of satellites when they move on near-polar orbits. 



7 SUMMARY AND CONCLUSIONS 

We coupled together two successful models for dynamical 
friction and tidal stripping and compared their predictions 
with high resolution N-body simulations to address the evo- 
lution and ultimate fate of satellite haloes within cosmic 
structures. Under the assumption that haloes are well de- 
scribed by an NFW profile we are able to predict if a satel- 
lite of given mass, orbital eccentricity and infall redshift, will 
merge, evaporate or survive under the simultaneous action 
of dynamical friction and tidal mass loss. 

We emphasize that we have obtained a complete pre- 
dictive scheme for the fate of a satellites whose masses at 
the time of infall into the main halo are known (below we 
refer to typical cosmological orbits): 

• High mass satellites (M St o > O.lMh) sink rapidly to- 
ward the center of the main halo without significant mass 
dispersal: The dynamical friction timescale for a rigid satel- 
lite (equation [25]) gives the correct timescale of merging. 

• For satellites of mass 0.01M h < M s ,o < O.lMh dynam- 
ical friction is still strong and drives the satellite toward 
the center where tidal mass loss become severe. Low con- 
centration satellites are disrupted, while high concentration 
satellites, severely pruned by the tidal field, survive with 
masses 0.01M Sj o, and settle into inner orbits with a typical 
reduction of the apocenter radii of a factor ~ 0.1 lower rel- 
ative to the initial value. The dynamical friction timescale 
for these satellites is longer than for their rigid counterpart, 
and is given by equation (27). 

• Light satellites with mass M s ,o < 0.01M h are almost 
unaffected by dynamical friction which is operating on a 
rather long timescale. Mass loss by the tidal field, which is 
not severe on these cosmological orbits, stabilizes further the 
orbit. 

• Low concentration satellites below O.lMh can be dis- 
rupted by tides before their orbital decay is complete. Com- 
parison of the dynamical friction timescale and the disrup- 
tion timescale as provided in this paper allows to find the 
actual lifetime of satellite haloes. 

We predict that, because of the combined action of 
stripping and dynamical friction, a primary halo at z = 
will host preferentially satellites with mass M s /Mh <C 0.01 
as the heavier ones would have been accreted or/and dis- 
persed in the background, leaving a "depression" in the mass 
function of substructure above 0.01 M s /Mh,(of course we are 
neglecting effects due to the evolution of the main halo it- 
self). This feature should be more evident in Milky Way-size 
haloes than in cluster haloes as in the former bound satel- 
lites had more time to evolve. 

Since the destructive power of the tidal field (and in par- 
ticular of tidal shocks) depends sensitively on the degree of 
circularity of the satellite orbit, a large galaxy halo like that 
of the Milky Way (> 1O 12 M ) should host satellites moving 
preferentially on circular orbits as a consequence of the selec- 
tive action of the tidal field. Also, because dynamical friction 
seems unable to render the satellites' orbit circular (van den 
Bosch et al. 1999; CMG99), the low eccentricities should 
have been already present as initial conditions. This "selec- 
tion effect" will be extremely weak for smaller satellites (be- 
low 0.01 — O.OOIMq) because their orbit barely decays and 
thus will have in general long survival times (only low con- 
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centration satellites could disappear quickly but they are not 
common in CDM models; see e.g. Eke, Navarro & Steinmetz 
2001 and Bullock et al.2001). This mass regime corresponds 
to that of the dwarf spheroidal satellites of the Milky Way. 
On the other end, the Magellanic Clouds, the dwarf ellipti- 
cal satellites of M31 and perhaps the dwarf spheroidal For- 
nax are all massive enough to fit in the intermediate regime 
where destruction is still possible; thus these galaxies could 
have survived because their host haloes had nearly circular 
orbits. In the case of the Magellanic Clouds a nearly cir- 
cular orbit is indeed measured (Kroupa & Bastian 1999). 
There is, however, at least one caveat to this interpretation, 
namely that both the dwarf ellipticals of M31 and the Mag- 
ellanic Clouds could be dense enough to survive shocks on 
even very eccentric orbits (Mayer et al. 2001b). Only when 
all the orbits of the satellites will be accurately determined 
we will know whether eccentricity or internal structure was 
more important in determining their survival. 

The calculations described in this paper can become a 
useful tool when coupled to cosmological simulations. The 
final goal is to find an appropriate description of the dynami- 
cal evolution of substructure in a halo. Increasing computing 
power and code performances has recently allowed to carry 
out extremely high resolution simulations that follow the 
evolution of substructure in dark matter haloes (Ghigna et 
al. 1998, 2000; Mayer et al. in preparation). These repre- 
sent the new ground where CDM models are being tested 
and their predictions compared to observations. However, 
these simulations remain costly and usually only one sys- 
tem at a time can be simulated down to very small scales. 
On the other end, resolving the mass function of substruc- 
ture in depth is important in light of the problem of the 
overabundance of satellites (Moore et al. 1999; Klypin et 
al. 1999). Such mass function can be viewed as the convo- 
lution of the mass function of satellites at an earlier epoch 
with an evolutionary filter function that depends on the dy- 
namical mechanisms analyzed in this work. Therefore, our 
results can allow to address the substructure problem in a 
statistical way orders of magnitude faster than with N-Body 
simulations; as an example one can explore a large number 
of dynamical histories by randomly varying the orbital and 
structural parameters in the range typical of cold dark mat- 
ter cosmogonies, and work is in progress (Taffoni et al. 2002 
in prep). Here we make a first attempt starting with uni- 
form distributions. Clearly the time dependent potential of 
the growing primary halo, whether it is a galaxy or a clus- 
ter, is an additional ingredient that only simulations can 
incorporate and which could affect the orbital dynamics of 
the satellites. However, the latter limitation can be partially 
overcome by using the merger tree extracted from a low 
resolution simulation, as done within some semi-analytical 
models (Somerville & Primack 1999; Kauffmann et al. 1999; 
Cole et al. 2000) or using analytical merger trees provid- 
ing a good approximation to the latter (Taffoni et al. 2002; 
Monaco et al. 2002). 

A key result of our analysis, and one that is in agree- 
ment with the high resolution cosmological simulations from 
which the initial orbits were drawn, is that the inner, most 
bound part of small satellites as concentrated as expected 
in CDM models (Eke, Navarro & Steinmetz 2001; Bullock 
et al. 2001) survive for timescales comparable to or longer 
than the age of the Universe. This residual has a size cor- 



responding to a few percents of the initial virial radius; this 
is comparable to the scale of the baryonic component in 
galaxies, so we can argue that galaxies will mostly survive 
within the main halo. This result is also confirmed by high- 
resolution SPH simulations of the formation of Milky Way- 
like galaxy (Governato et al. 2002). Indeed, dissipation could 
make the inner part of the haloes even more robust against 
tides (Navarro & Steinmetz 2000). On the other end, ad- 
ditional tidal shocks occurring during encounters between 
substructure haloes, i.e. galaxy harassment (Moore et al. 
1996, 1998), might have a counteracting effect and could ac- 
tually increase mass loss. However, detailed simulations of 
this mechanism have shown that only very fragile, LSB-type 
galaxies would be severely damaged by harassment (Moore 
et al. 1999); halo profiles of these galaxies likely correspond 
to the low concentration satellites studied in this paper (Van 
den Bosch & S waters 2001) which we have shown are easily 
disrupted even by the tides of the primary halo alone. Thus, 
adding harassment would only accelerate the disruption of a 
few satellites while not affecting the survival of the majority 
of them which, in CDM models, have high concentrations. 
Hence the picture emerging from the life diagrams of the 
satellites shown in this paper is robust. Satellites close to 
disruption at the present time, like Sagittarius in the Milky 
Way subgroup, must have been much bigger in the past in 
order for dynamical friction to drag them to an inner orbit 
where dissolution can easily take place; alternatively they 
could have entered the halo at fairly high redshift, which 
would place them naturally on an inner, tightly bound orbit 
(Mayer et al. 2001b). In clusters, dwarf galaxies cannibalized 
by giant CD galaxies might also trace an early population. 

Satellites infalling at redshift one or lower in the 
main haloes will complete several orbits and eventually un- 
dergo morphological changes by tidal stirring (Mayer et al. 
2001a,b) and harassment (Moore et al. 1996, 1998). These 
will produce diffuse streams of stars while they are orbiting 
(Helmi & White 1999; Johnston, Sigurdsson & Hcrnquist 
1999), contributing to the build up of an extended stellar 
halo population. Such population should be present out to 
more than 200 kpc in the Milky Way halo, as the plung- 
ing orbits of satellites seen in cosmological simulations go 
this far out. On the contrary, a less extended stellar halo 
should be expected if dynamical friction were more efficient 
in dragging satellites to the center. The amount of stellar 
halo substructure out to large distances could thus reveal the 
original mass function of observed dwarf spheroidal galaxies 
in the Local Group. Components decoupled in their kine- 
matics as well as in the metallicity and age of their stars 
should be present, but tracking such properties might be a 
daunting task observationally if enough phase mixing occurs 
(Ibata et al. 2001a,b); however, while in the inner halo fast 
orbital precession and heating by other clumps might blur 
the streams, the phase space distribution of the outer halo 
material should still carry the memory of the initial orbits 
of the satellites (Mayer et al. 2002). 

The prescriptions for the decay and disruption rates 
obtained in this work provide a complete framework which 
can improve the predictive power of semi-analytical models 
of galaxy formation. They enable to follow the complex evo- 
lution of substructures in hierarchical models in a straight- 
forward manner. 
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APPENDIX A: CALCULATION OF THE TIDAL 
ENERGY FOR A NFW PROFILE 

At each pericenter passage the satellite crosses very rapidly 
the central and more concentrated regions of the primary 
halo. The duration of those encounters is fast compared with 
the dynamical time of the object. Such kind of interactions 
are called tidal shocks (Spitzer 1987). We will use the results 
derived by GHO to describe the amount of heating due to 
tidal shocks on a satellite moving inside an extended mass 
distribution. 

During an orbital period P or b the tidal force f s , t id per 
unit mass produces a global variation on the velocity of the 
internal fluid: 

f Por b 

Av= / t tid dt, (Al) 
Jo 

where we have applied the impulse approximation in the hy- 
pothesis that the time scale of interaction is short compared 
with the dynamical time of the satellite (t = refers to the 
initial satellite position at apocenter distance). 
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In a spherically symmetric system of mass Mh, the tidal 
force per unit mass exerted by the background on a dark 
matter particle of the satellite is: 



ftid = 



4(3/1 - //)(? ■ R s )f - jiiRs] , 



(A2) 



where r = r/i?h is the direction to the center of mass of the 
satellite (CMS), R s is the position of the particle respect to 
CMS. Note that is the virial radius of the main system. 
Here: 

M(r) 



fj,(r) 



is the adimensional mass profile, and: 
dfi(r) 



n'(r) 



d In r 



(A3) 



(A4) 



For a NFW profile [i and fi 1 axe functions of the normalized 
radius x = r/Rh and of the concentration Ch of the primary 
halo: 

fi(x, Ch) = 

and 



ln(l + c h x) — c h x/(l + c h x) 
ln(l + c h ) - c h /(l + c h ) 



ln(l + c h ) 



1 / c h x y 

- Ch/(1 + Ch) VI + ChXj 



(A5) 



(A6) 



In the case of stable orbits the angular momentum J is 
conserved and we can use the identity 

dt = (r 2 /J)dO 



to re- write equation(Al) into components (GHO): 
Av = -^{{Bi - B a )x, (B 2 - B a )y, -B 3 z} , 



where 

Si(ch) = 

B 2 (c h ) = 

B 3 (c h ) = 



Fi(x, Ch) cos 9d8 

Fi(x,Cb) sin 2 6d6 
M^£h) 



(A7) 



(A8) 



(A9) 



(A10) 



(All) 



where 9 m is the maximum value of the position angle, and: 
Fx{x,Cb) = 3 

[ln(l + Chx) - c h x/(l + c h x)] - [c h a:/(l + c h x)] 2 
X a;[ln(l + c h ) - c h /(l + Ch)] ' 1 ' 

This velocity changes cause a reduction of the binding en- 
ergy of the system: 



(AE) = ( -|Av| 



(A13) 



Averaging over an ensemble of dark matter particles in 
a spherically symmetric satellite we have that (a: 2 ) = (y 2 ) = 
(z 2 ) = R 2 I 2, and using equation(A8), the tidal energy 
gained by the satellite becomes: 



<-> - (ft)" 



(B 1 - B 3 ) 2 + (B 2 - B 3 ) 2 + B 2 



(A14) 



1.5 - 



0.5 




Figure Al. The intensity of the tidal force |ftid(t)| normalised 
to its value at the first periastron |ftjd perl- We plot the module of 
the tidal force (equation [A2]) as a function of time for a satellite 
of mass M s , = 0.01Af h and e = 0.7 x c (E) = 0.5. The solid 
line refers to a stable orbit, the dashed line to an unstable one. 
When the drag force is active, the intensity of the tidal force, and 
consequently of the shock energy, grows with time. 



In the previous expression, the contribution due to the halo 
and the orbital parameters (e and x 2 [E\) is confined in the 
function: 

-r tv\ i { GM h V 1 

^[c h ,x c (E),e] = 



(Bi - B 3 ) 2 + (Ba - B 3 ) 2 + Bl 



(i 



A15) 



here Vh is the circular velocity of the main halo at virial 
radius. 

It is then useful to write the shock energy as: 

(AE) =S[c h ,x c (E),e] R 2 . (A16) 

When the frictional drag force is active, it is not possible to 
change the integration variable according to equation (A7). 
The energy change becomes: 

<A*> - m 



(A, - A 3 ) 2 + (A 2 - A :i ) 2 + A% 



Here 

Ai{e,x c [E]) = 

A 2 (e,x c [E]) = 

A 3 (e,x c [E]) = 
with 



F 2 (x, c h ) cos 9 dt 



F 2 (x, Ch) sin 2 9 dt 



A*(x,Ch) 



dt, 



R 2 s ■ (A17) 

(A18) 
(A19) 
(A20) 
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F 2 (x,c h ) 



Fi(x, Ch) 



(A21) 



Once again we separate the contribution due to the environ- 
ment: 

2 

i i : \ i, \ 

F[c^x c {E),e] 



(GMhX 

{ Rl ) 

(Ai - A 3 ) 2 



(A 2 - A 3 ) 2 + Al 



;A22) 



For an unstable orbit 



{AE) = f[c h ,x c (E),e] i? 2 



(A23) 



The shock energy in this case must be evaluated along the 
perturbed orbit. As the drag force drives the satellite in the 
internal region of the halo the (AE) increases (Fig. [Al]). 



tdis 



N 



(B8) 



This formula provides a simple estimate of the disruption 
time valid on cosmological relevant orbits with a precision 
of 25 per cent. 



APPENDIX B: AN APPROXIMATE ESTIMATE 
OF THE DISRUPTION TIMESCALE 

Because the lifetime of light satellites is mostly set by tidal 
disruption we estimate here the disruption time. If the main 
halo density profile is isothermal (ISO) then GHO showed 
that the shock energy change is: 

/ap\ ^ M 2 2sin 2 6> m + 4fl 2 n 2 

{AE)lS ° = UJ 6(exe) 2 A{Xt] ■ R * ' (B1) 

where # m is the maximum value of the position angle which 
varies from tv/2 to n. Using the orbit equation (equation [6]) 
we can evaluate # m : 



r apo ^ 

= 2ex c - (B2) 

"Der 



R f r 



dx 

x 2 y / ln(r c /x) 2 — (r c /x) 2 e 
and the orbital period 

^ ap ° _dx 

yj\n(r c /x) 2 — (r c /x) 2 e 2 — 1 



The shock in the ISO profile equals the shock of an 
NFW case when Ch = 30. We have than: 

(A£)nfw ~ (AE) lso x (0.029 c h + 0.13) . (B4) 

At each pericenter passage the satellite is shock heated 
and its radius i? s is reduced of a factor AT?. As an approxi- 
mation A_R ~ Rs,o/N where iV is the number of pericenter 
passages necessary to destroy the satellite. Then, we have 
an implicit equation for iV 

JV-l 

I\ * — ' (Afi)NFW,0 
i=l 

where E — 0.5GM B /R By0 , and (A_E) N fw,o is evaluated at 
the initial half mass radius. Since 



JV-l 



£Y = N(N-1H2N-1) (B6) 

i=l 

for large N equation (B5) would give: 



N~J48-—^ . (B7) 

The disruption time can then be written as: 
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